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ABSTRACT 

We present a determination of the distributions of photon spectral index and gamma-ray flux - the 
so called LogA^-Log6' relation - for the 352 blazars detected with a greater than approximately seven 
sigma detection threshold and located above ±20° Galactic latitude by the Large Area Telescope of 
the Fermi Gamma-ray Space Telescope in its first year catalog. Because the flux detection threshold 
depends on the photon index, the observed raw distributions do not provide the true hogN-LogS 
counts or the true distribution of the photon index. We use the non-parametric methods developed by 
Efron and Petrosian to reconstruct the intrinsic distributions from the observed ones which account 
for the data truncations introduced by observational bias and includes the effects of the possible 
correlation between the two variables. We demonstrate the robustness of our procedures using a 
simulated data set of blazars and then apply these to the real data and find that for the population 
as a whole the intrinsic flux distribution can be represented by a broken power law with high and low 
indexes of -2.37±0.13 and -1.70±0.26, respectively, and the intrinsic photon index distribution can be 
represented by a Gaussian with mean of 2.41±0.13 and width of 0.25±0.03. We also find the intrinsic 
distributions for the sub-populations of BL Lac and FSRQs type blazars separately. We then calculate 
the contribution of Fermi blazars to the diffuse extragalactic gamma-ray background radiation. Under 
the assumption that the flux distribution of blazars continues to arbitrarily low fluxes, we calculate 
the best fit contribution of all blazars to the total extragalactic gamma-ray output to be 60%, with a 
large uncertainty. 

Subject headings: methods: data analysis - galaxies: active - galaxies: jets - BL Lacertae objects: 
general 



1. INTRODUCTION 

A vast majority of the extragalactic objects observed 
by the Large Area Telescope (LAT) on the Fermi 
Gam ma-ray Space Telesc ope can be classified as blazars 
(e.g. lAbdo et all l2010a[ ). a unique subclass of active 
galactic nuclei (AGNs) for which t he jet is aligned with 
the o bserver's line of sight (e.g. iBlandford fc Konigll 
1 1 9 79h . Analyses of the gamma-ray spectra of blazars 
along with other signatures of AGNs indicate that the 
gamma-ray emission is an essential observational tool 
for understanding of the physics of the central engines 
of AGNs. In addition, as in all AGNs, the distribution 
of spectral and other characteristics of blazars, and the 
correlations among these characteristics and their cosmo- 
logical evolutions, are essential information for the stud- 
ies of the for mation and gro wth of central black holes of 
galaxies (e.g. llDermer Il2007h . 

This information comes from the investigation of the 
population as a whole. The process for any extragalactic 
source starts with the determination of the LogTV— Log^ 
relation which can be carried out simply by counting 
sources even before any redshifts are measured and dis- 
tances are used to determine intrinsic characteristics such 
as luminosities and source densities (and their evolution). 
Although redshifts are measured for many blazars the ex- 
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tant sample is not yet sufficiently large to allow an accu- 
rate determination of the intrinsic characteristics. Our 
ultimate goal is to carry out such an analysis but the 
focus of this paper is the determination of the flux and 
photon spectral index distributions. 

The detection threshold flux of blazars by the Fermi- 
LAT depends strongly on an object's gamma-ray spec- 
trum, such that harder spectra are detected at lower 
fluxes (measured for a given photon energy, here for pho- 
tons > 100 MeV). This means that for determination of 
the flux distribution we need both a measure of flux and 
the photon index T, and that one deals with a bi-variate 
distribution of fluxes and indexes, which is truncated be- 
cause of the above mentioned observational bias. Thus a 
bias free determination of the distributions is more com- 
plicated than just counting sources. 

The re have been analyses of this data (e.g. lAbdo et al.l 
2010b) using Monte Carlo simulations to account for the 
detection biases. In this paper we use non-parametric 
methods to determine the distributions directly from the 
data at hand. As stressed by Petrosian (1992), when 
dealing with a bi-(or more generally multi)-variate distri- 
bution, the first required step is the determination of the 
correlation (or statistical dependence) between the vari- 
ables, which cannot be done by simple procedures when 
the data is truncated. We use the techniques developed 
by E fron and Petrosian (EP, lEfron fc Petrosianl |l992, 
1999) which can account reliably for the complex obser- 
vational selection biases to determine first the intrinsic 
correlations (if any) between the variables and then the 
mono-variate distribution of each variable. These tech- 
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Fig. 1. — Flux and photon spectral index for the 352 Fermi- 
LAT blazars used in this analysis, those with test statistic > 50 
and |b| > 20°. BL Lac type blazars (n=163) are shown as blue tri- 
angles, FSRQs type blazars (n=161) are shown as red plus signs, 
and blazars of unidentified or ambiguous type (n=28) are repre- 
sented by black x's. It is seen that there is a selection bias against 
soft spectrum sources at fluxes below ~ 10 -7 photons cm -2 sec -1 . 
We also show for a selection of sources (but only a few for clarity) 
the approximate limiting flux for that source - that is the lowest 
flux it could have and still be sufficiently bright to be included in 
the sample given its location on the sky given the reported detec- 
tion significanc e. T he location of the line used for the truncation 
boundary (see 33. 1[) is shown in Figure [8] 

niques have been proven useful for application to many 
sources with varied characteristics and mos t recently to 
radio and optical luminosity in quasars in iSingal et al.l 
(|2011h . where a more thorough discussion and references 
to earlier works are presented. 

In this paper we apply these methods to determine 
the correlation and the intrinsic distributions of flux and 
photon index of Fermi-L AT blazars. In ^21 we discuss the 
data used from the LAT extragalactic catalog, and in $3] 
we explain the techniques used and present the results. 
In $4]we describe how the result from such studies are im- 
portant for understanding the origin of the extragalactic 
gamma-ray background (EGB) radiation. A brief dis- 
cussion and summary is presented in §5. A test of the 
procedures using simulated data set is discussed in the 
Appendix. 

2. DATA 

In this analysis we use the sources reported in the 
Fermi- LAT first y ear extragalactic source catalog (e.g 
lAbdo et al.l l2Q10ch . In particular we rely on the subset 
of sources that have a detection test statistic TS > 50 
and which lie at Galactic latitude \b\ > 20°. This is the 
same criterion adopted by the LAT t eam fo r analysis of 
the blazar population (jAbdo et al.l (j2Q10bh - hereafter 
MA) and includes those sources that are fully calibrated 
and removes spurious sources. The test statistic is de- 
fined as TS = — 2 x (ln(Lo) — ln(Li)), where Lq and 
Li are the likelihoods of the background (null hypoth- 
esis) and the hypothesis being tested (e.g., source plus 
background). The significance of a detection is approx- 
imately n x a = y/TS . Of 425 total such sources, 352 
are identified as blazars. The rest are either identified as 
radio galaxies (2), other AGN or starbursts (6), high lat- 
itude pulsars (9), and objects without radio associations 
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Fig. 2. — Correlation factor (3 versus test statistic r for a photon 
index and flux correlation of the form given in Equation [2] for the 
352 blazars used in this analysis (solid curve), the subset of BL Lac 
type blazars (dashed curve), and the subset of FSRQs type blazars 
(dash-dot curve). The la range of best fit values for (3 are where 
|r| < 1. For comparison, the dotted curve shows the correlation 
factor for just those sources above 5 x 10 -8 photons cm -2 sec -1 , 
where the data truncation in the Fioo,r plane is not as relevant. 

(56). Among the blazars 161 are identified as Flat Spec- 
trum Radio Quasar (FSRQs) type, 163 are identified as 
BL Lacertae (BL Lac) type, and 28 have uncertain type. 
The fluxes and photon indexes of the blazars are plotted 
in Figure [TJ 

The 352 blazars used in this analysis range in gamma- 
ray flux (integrated over the photon energy range 100 
MeV to 100 GeV from a power law fit to the Fermi- 
LAT data and designated here as Fioo) from 9.36 x 10 -9 
to 1.37 x 10 -6 photons cm -2 sec -1 . The photon index 
T, obtained by fitting a simple power-law to the spec- 
tra in the above energy interval, ranges from 1.253 to 
3.039. The photon index T is defined such that for the 
monochromatic photon spectral density n(E)dE oc E~ T 
(or the vF v oc v~ T+2 ). The bias mentioned above is 
clearly evident; there is a strong selection against soft 
spectrum sources at fluxes below Fioo ~ 10 -7 photons 
cm -2 sec -1 , caused by the dependence of t he Fermi-LAT 
point spread function (PSF) with energy (|Atwood et al.l 
l2009h . 

Each source has a TS associated as discussed above, 
and the backgroun d flux is a function o f position on the 
sky, as discussed in lAbdo et"aT1 (|2010d ) . In Figure Q] we 
also show the approximate limiting flux of some (not all 
to avoid confusion) objects, an estimate of the lowest flux 
it should have (at its location in the sky and having the 
specific value of its index) to be included in the sample, 
given by Fu m = Fioo/ y/TS/50. However, as discussed 
below, because the limiting flux as determined in this way 
is not the optimal estimate, we use a more conservative 
truncation as shown by the straight line in Figure [51 

3. DETERMINATIONS OF DISTRIBUTIONS 

3.1. Correlations 

As stressed above, when dealing with a bi-variate 
truncated data it is imperative to determine whether 
the variables are independent or not. If flux and pho- 
ton index are independent, the combined distribution 
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G(Fioo,r) can be separated into two independent distri- 
butions V{^ioo) and h{T). However, independence may 
not be the case for Fi o and T even though flux is a 
distance dependent measure while photon index is not. 
An intrinsic correlation between the photon index and 
luminosity may be strong enough to manifest as a flux- 
index correlation even after cosmological smearing of the 
correlatio ns, as is the case for example in ga mma-ray 
bursts (jYonetoku et alJl2004t ILlovd et al.ll2000t ). Even if 
there is no intrinsic correlation between the photon index 
and flux, if the selection process introduces some corre- 
lation then the independence assumption breaks down, 
and any such correlation should be removed in order to 
obtain bias free distributions of the variables Thus, the 
first task is to establish whether the variables are inde- 
pendent. Determining the correlation when the data is 
truncated is not straight forward. 

We use the Efron-Petrosian method to determine 
whether the two variables are correlated. This method 
is a version of the Kendall Tau statistic test devised for 
truncated data and uses the test statistic 

to test the independence of two variables in a data set, 
say (xj,jjj) for j = 1, . . . , n. For untruncated data (i.e. 
data truncated parallel to the axes) IZj is the y rank of 
the data point j within the set with Xi < xj (or alter- 
natively Xi > Xj), which we call the associated set. If 
the data is truncated, say it includes only points with 
V > 2/iim = d( x ) then the associated set is defined as the 
largest untruncated set of points associated with x^ i.e. 
not all points Xi > Xj but only a subset of these that 
have yk > yiimj = tfoj) ( see EP for a full discussion of 
this method). 

If (xj,yj) were independent then the rank IZj should 
be distributed uniformly between and 1 with the ex- 
pectation value and variance Sj = (1/2) (j + 1) and 
Vj = (1/12) (j 2 + 1), respectively. Independence is re- 
jected at the n a level if | r | > n, if r turns out to be 
significantly different than expected value of zero. In 
such a cast the correlation is removed parameterically as 
follows. We define a new variable y r = f(x,y) and repeat 
the rank test with different values of parameters of the 
function / and determine the nature of the correlation 
by the best fit value of the parameters that give r = 0; 
the na range is obtained from — n < r < n. 

We carry out this test for our data set using a variable 
transformation, which is a simple coordinate rotation, by 
defining a new variable we call the "correlation reduced 
photon index" as 

r cr = r - a x log (^2° ) . (2) 

and determine the value of the parameter j3 empirically 
that makes Fioo and T cr independent, which then means 
that the distributions of Fioo and T cr are indeed separa- 
ble: 

4 The observed distribution (in Figure [T} clearly shows a strong 
correlation. However, most of this correlation is due to the data 
truncation described above. 



G(F 10 o,r) = v<F 10 o) x Hr CI ). (3) 

Once the monovariate distributions are determined then 
the true distribution of T can be recovered by an inte- 
gration over Fioo as: 

m = 

jf HF 100 )h \T-f3x log (^)) dF W0 . (4) 

Here Fq is some fiducial flux we chose to be Fo = 6 x 10 -8 
photons s _1 cm -2 sr _1 which is approximately where the 
flux distribution breaks (see below), although its value is 
not important. The data described in £j2]are truncated in 
the Fioo — T plane, due to the bias against low flux, soft 
spectrum sources. We can use a curve approximating 
the truncation, Tii m = g(Fioo), which allows us to define 
the associated sets. The associated set for each point are 
those objects whose photon index is less than the limiting 
photon index of the object in question with its specific 
value of Fioo- 

We have tested this procedure using a simulated 
dataset from the Fermi-LAT collaboration designed to 
resemble the observations, but with known distributions 
of uncorrelated photon index and flux and subjected to a 
truncation similar to the actual data. The results are de- 
scribed in the Appendix where we demonstrate that we 
can recover the input distributions which are of course 
quite different than the observed biased distributions. As 
shown in the Appendix we find that we recover the input 
distribution best if we start with a truncation boundary 
Thm = g(logFioo) roughly defined by the limiting values 
obtained from the TS values. We then gradually move 
this limit to higher fluxes (see Figured]) and to more con- 
servative estimations of the truncation. This procedure 
is stopped when the results do not change significantly. 
This way we lose some data points but make certain that 
we are dealing with a complete sample with a well de- 
fined truncation. Note that when defining new variables 
the truncation curve as a function of flux should also be 
transformed by same parameter /3; 

r c r,iim = r lim -p x Log ( - 100 ) . (5) 

\ri00-min J 

We subject the actual data to the same procedure, and 
the results converge with the same cutoff limit location. 
Figure [2] shows the result of the test statistic r as a func- 
tion of the correlation parameter j3 for a all blazars and 
the subsets including only BL Lacs and FSRQs in the 
sample. Table 1 shows the best fit values and 1<j ranges 
of the correlation parameter j3. We note that the corre- 
lation is weak, for example for all Fermi blazars the best 
fit value is /3=0.02 ± 0.08 indicating a weak correlation 
with the la range including f3 = or no correlation. 
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Fig. 3. — Observed (diamonds) and intrinsic (stars) cumula- 
tive distribution of flux <E>(Fioo) = f£S ^^100)^100 f° r tne 352 
Fermi-h AT blazars used in this analysis, shown for the best fit 
value of the correlation parameter (3. The error bars represent the 
la range of the correlation parameter (3, and are in general smaller 
than the stars, and are larger tha n th e statistical error. The nor- 
malization is obtained by equation 1101 
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Fig. 4. — Observed (diamonds) and reconstructed intrinsic 
(stars) differential distribution of flux V{^ioo) for the 352 Fermi- 
LAT blazars used in this analysis. The error bars represent the 
la range of the correlation parameter (3. The intrinsic distribu- 
tion is a power law with a break at F^ r ~ 6 x 10 -8 photons 
cm -2 sec -1 . The best fit slopes for the intrinsic distribution are 
-2.37±0.13 above the break and -1.70±0.26 below, and the best fit 
intrinsic distribution is plotted as the dotted line. We also plot 
V{^ioo) as determined in MA (small crosses), with error bars (dot- 
ted lines). The best fit value for ^F break ) is 2.2 xlO 8 sr -1 ^^. 

3.2. Distributions 
With the correlation removed the independent distri- 
butions V(^ioo) and ft(T cr ) can be d etermined using a 
method outl i ned in iPetrosianl (jl992h and developed by 
iLynden-Belll (jl971[ ). These methods give the cumula- 
tive distributions by summing the contribution from each 
point without binning the data. 

3.2.1. flux distributions 
For the flux the cumulative distribution 
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Fig. 5. — Intrinsic cumulative distribution of photon index 
t\T cv ) = Jq CT h(T CY )dT CY for the 352 Fermi-LAT blazars used in 
this analysis. The normalization of P is arbitrary 
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Fig. 6. — Observed (diamonds) and reconstructed intrinsic (stars) 
distribution of photon index h(T) for the 352 Fermi-LAT blazars 
used in this analysis. The intrinsic distribution is calculated from 
the flux distribution and the correlation reduced photon index dis- 
tribution by equation 0] The stars represent the intrinsic distribu- 
tion calculated with the best fit value of the correlation parameter /3 
and the solid curve is the best fit Gaussian function to these values, 
while the dotted curves represent the best fit Gaussian functions to 
the extremal intrinsic distributions allowed by the la range of (3. 
The intrinsic distribution can be represented by a Gaussian with 
a mean of 2.41±0.13 and la width of 0.25±0.03, while the ob- 
served distribution can be represented by a Gaussian with a mean 
of 2.32±0.01 and la width of 0.32±0.01. The normalization of h(T) 
is arbitrary. 



is obtained as 



^100) =n i 1 



(7) 



where j runs over all objects with fluxes -Fioo,j > -^100? 
and N(J) is the number of objects in the associated set 
of object j; namely those with with a value of Fioo,i > 
Fiooj and r cr) i < r crj ii m (Fioo,j) determined from the 
truncation curve described above. Equation [71 represents 
the established Lynden-Bell method. We note again that 
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the cutoff curve as a function of flux is scaled by /? in the 
same manner of equation O The use of only the asso- 
ciated set for each object removes the biases introduced 
by the truncation. 
The differential distribution 



^100) = 



gLFioo 



(8) 



is obtained by fitting piecewise polynomial functions via 
least-squares fitting to <£(Fioo) and calculating its deriva- 
tive. Figure [3] shows the true intrinsic and (the raw) 
observed cumulative distributions of Fi o for all 352 
blazars, while figure |4] shows the calculated true intrin- 
sic differential distribution VO'ioo)? along with those ob- 
tained from the raw observed data without correcting for 
the bias. A direct comparison to the results from MA is 
presented there as well. The differential counts manifest 
a broken power law which can be fit by the form 



^100 ) = HFbreak) 



"100 



break ) 



F breai 



for Fioo > F bre $) 



Fi 



break 



for Fioo < F h r 



ik • 



mabove and rribeiow are the power law slopes above and 
below the break, respectively, and are obtained from a 
least-squares fitting of t/(Fioo), 

as is the value of Fbreak • 
At values of Fioo above F^t = 1 x 10 -7 photons cm -2 
sec -1 the truncation is not significant and we can obtain 
the normalization by scaling the cumulative distribution 
$(^100) such that. 



$(F NT ) 



n±Vn 

8.26 sr ' 



(10) 



where N is the number of objects above F^t and is equal 
to 60 for all blazars, 12 for BL Lacs, and 48 for FSRQs. 
The \fN uncertainty arises because of Poisson noise for 
the brightest sources, and 8.26 sr is the total sky coverage 
considered, which is \b\ > 20° as discussed in $21 This 
also gives the value of i^Fbreak) by inegrating V0*ioo) at 
fluxes above Fnt and setting this equal to $(Fnt)' 



break ) 



$(Fatt) {m above - 1) F N ™ above 1 F^ ak 



(11) 

The best fit value for ^i^rea/e), corresponding to the 
best-fit value of m a bove, is then 2.2 xlO 8 sr -1 ^ - ^. 
3.2.2. Photon index distributions 

A parallel procedure can be used to determine the 
distribution of the correlation reduced photon index, 
namely the cumulative distribution 



kr' cr )dr' CI 



obtained with 



^r cr ) = n(i 



M(k) 



(12) 



(13) 



in the method of lLvnden-Belll (|1971[ ) can be differentiated 
to give the differential distribution 



HTcr) = 



dP(T cr ) 
dT cr 



(14) 



In this case, k runs over all objects with a value of 
r cr ,k < T cr , and M(k) is the number of objects in the 
associated set of object fc; i.e. those with Y cr ^ < T cry k 
and Fioom > Fu m ^ obtained from the truncation line at 

rcr,k- 

Figure [5] shows the cumulative distribution of the cor- 
relation reduced photon index P(r cr ) for all 352 blazars. 

Differentiation of this gives /i(r cr ), which can be substi- 
tuted in equation |4] to obtain the intrinsic distribution 
of the photon index itself, h(T). The results are shown 
in Figure [6] along with the raw observed distribution for 
comparison. Because the mean of intrinsic distribution 
of photon index is sensitive to the value of the correlation 
parameter /3, we include the full range of intrinsic distri- 
butions resulting from the la range of /?. A Gaussian 
form provides a good description of the intrinsic distri- 
bution of the index. 

We have carried out identical procedures to obtain the 
distributions of the BL Lac and FSRQ subsets of the 
data. Table 1 summarizes the best fit parameters for 
the intrinsic flux and photon index distributions, for the 
sample considered as a whole, and for the BL Lac and 
FSRQs sub-populations separately. The errors reported 
include statistical uncertainties in the fits and the de- 
viations resulting from the la range of the correlation 
parameter j3. A higher value of j3 (i.e. more positive cor- 
relation between flux and photon index absolute value) 
moves the mean of the photon index distribution down 
to a lower absolute value of the photon index and makes 
the faint end source counts slope less steep (less negative 
rribeiow), while a lower value of j3 has the opposite effect. 

3.3. Error Analysis 

It addition to uncertainty in the value of j3 and those 
due to the fitting procedure there are other effects that 
can add to the uncertainties of the final results. Here 
we consider the effects of some factors which we have 
ignored in the above analysis. 

1. Individual measurement uncertainties: We have 
treated individual sources as points having a delta func- 
tion distribution in the fl ux-index plane, r esulting in a 
possible Eddington bias (Ed dingtonl 119401 ). The mea- 
surement uncertainties can be included by changing the 
delta functions to kernels whose widths are determined 
by the reported measurement errors. The main effect of 
this will be smearing out of the distribution which can be 
neglected if the errors are small compared to the width 
of features in the distributions. This effect will not in- 
troduce any bias as long as measurement uncertainties 
are symmetrical about the reported value (e.g the ker- 
nels are Gaussian) and the distributions themselves are 
symmetrical or fairly flat . The forme r is the case for the 
reported uncertainties in lAbdo et"al] (j2010af ) . This later 
is the case for the distribution of T, where the reported 
measurement errors vary between 0.04 and 0.35, but this 
is not likely to introduce any bias given the symmetrical 
distribution in V. 

The situation is different for the flux distribution, 
which is a power law. In this case the typical flux un- 
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certainty values are on the order of 1/4-1/3 of the re- 
ported fluxes. There may be more or fewer sources be 
included than missed in a flux limited sample such as 
this one. For example, for a power law differential distri- 
bution with index |m& e £ow| > 1, which is the case here, 
more sources will be included than missed at any flux 
for which there are errors in the reported fluxes, which 
will bias the distribution. The effect is different for the 
case where fractional measurement errors are constant 
with flux versus when they change with flux. For con- 
stant fractional flux measurement errors, an error will 
be introduced on the normalization of the source counts 
and can be appro ximated by [1/2 S 2 vn\) e i ow (mb e i ow + 1)] 
(jTeerikor pi 2004) where 5 is the fractional error in flux. 
For 

ir^beiow ~ — 1-7 and S ~ 0.3 this error will be about 
5%, which is small compared with other sources of nor- 
malization uncertainty. On the other hand, there will be 
an effect on the reconstructed slope of the counts only 
if the fractional flux measurement errors change with 
flux. It is expected that the fractional flux errors will 
be larger for lower fluxes, and in the extreme case that 
they do increase from negligible at high fluxes to values 
of 1/4-1/ 3 at the lowest fl uxes, according to simulation 
results in iTeerikorpil ([20041 ) resulting fractional errors in 
the source count slope will be around 7%, which is signif- 
icantly smaller than the errors we already quote for the 
s ource count slopes. 

iCaditz fc Petrosianl (j 19931 ) evaluated the effect of mea- 
surements error on EP method determinations of lumi- 
nosity functions of quasars using a Gaussian kernel and 
found that for individual uncertainty widths significantly 
smaller than the data range, the effects of the inclusion 
of measurement uncertainties were small (e.g. Figure 1 
of that work). Given the Gaussian symmetrical nature 
of the reported uncertainties, the symmetrical nature of 
the photon index distribution, and the relatively shallow 
faint end power law slope of the source counts distribu- 
tion, and especially the relative size of the reported un- 
certainties compared to the range of values considered, in 
light of the analysis done by Caditz and Petrosian we con- 
sider the errors introduced by the individual data point 
uncertainties to be negligible compared to the uncertain- 
ties introduced by the range of the correlation parameter 
j3 and the uncertainties in the power-law fits. 

2. Blazar variability: It is well known that blazars 
are inherently variable objects. There are two potential 
effects arising from blazar variability relevant to the anal- 
ysis here. One is similar to measurement error discussed 
above, in that it presumably would cause more objects 
to rise above the flux limit and be included in the survey 
than go below the flux limit and be excluded. The other 
is that the reported temporally averaged quantities such 
as flux and photon index, which we use in this analysis, 
may deviate from the true average values. 

Addressing the former issue, as discussed in the first 
year Ferm i-LAT extragalacti c source catalog (e.g. Fig- 
ure 11c of lAbdo et al.ll2QlQch . the pattern of maximum 
flux versus mean flux does deviate at the lowest detected 
fluxes. This indicates that variability becomes more im- 
portant with decreasing flux, but not as sharply as might 
be expected from previous EGRET data. This will be 
even less sharp for the TS > 50 sources we use as op- 
posed to the entire TS > 25 s ample considered there. As 
shown in lAbdo et al.l (|2010ef h the peak-to- mean flux ra- 



tio is a factor of two or less for most Fermi-LAT blazars, 
which excludes the possibility that most of the sources 
are detected because of a single outburst which happened 
during the 11 months of observation and are undetected 
for the remaining time. We believe the bias resulting 
from detecting blazars only in their flaring; state is small . 

A ddressing the late r issue , both lAbdo et al.l (|20Tot ) 
and Ack ermann et al.l (|2011f ) presented a detailed anal- 
ysis of the variability issues with Fermi-LAT blazars. 
They find that most sources exceed their average flux for 
less than 20%, and often less than 5%, of the monitored 
time, and conclude that both the timescale of variability 
is short compared with the length of observations, and 
that the measured average quantities are n ot highly bi- 
ased by flaring. Moreover, as also shown in lAbdo et al.l 
(|2010ef ). there is little or no temporal variation of the 
photon index with flux. We thus believe that no large 
systematic uncertainties result from the use of these av- 
eraged physical quantities 

3. Source confusion: Source confusion can also intro- 
duce errors at the faint end of the reconstructed distribu- 
tions because of relatively broad point spread function of 
the Fermi-LAT; some faint sources may be either missed 
entirely or erroneously combined into the fluxes of other 
sources. We first note that these two phenomena will 
have opposite systematic effects on the faint end source 
counts slope, as the former would tend to make it less 
steep while the later would tend to make it steeper. In 
addition to this self-canceling tendency, several tests ar- 
gue against Fermi-LAT' s blazar detec t ions be ing signif- 
icantly confusion limited. lAbdo et al.l (|2010al ) estimates 
that at Galactic latitudes above ±10° and at a TS > 25 
detection threshold, approximately 7.6% of sources (80 
out of 1043) are missed because of confusion, and blazars 
are 85% of the \b\ > 10° sources. Since we have used 
only sources detected at TS > 50 and at latitudes above 
±20° Galactic latitude, the effect of confusion should be 
lower because the sample will be more complete. Again, 
the faint end source counts slope could be altered by 
an amount considerably less than this. Other evidence 
against source confusion being a significant problem for 
Fermi-LAT blazars is the large increase in the number of 
extragalactic sources fro m the first-year to sec ond-year 
Fermi-LAT catalogs (e.g I Ackermann et al.l 120 ill ). Addi- 
tionally, there is the analysis described in section 8.3 of 
MA where different extragalactic sky scenarios were sim- 
ulated and run through an instrument detection and cat- 
alog pipeline, including an extremal scenario with a sin- 
gle steep power law distribution of blazars with a differ- 
ential slope of -2.23, in which blazars with fluxes greater 
than 10 -9 photons cm 2 s _1 would produce 70% of the 
total extragalactic gamma-ray background. Even under 
this much more dense sky scenario, many more blazars 
would be detected by the instrument and analysis, in- 
cluding at low fluxes. 

We also note that to the extent that the effects of both 
blazar variability and source confusion will have some 
photon index dependence, because of differing spectra 
in the case of the former and the Ferrm-LATs energy- 
dependent point spread function in the case of the later, 
then any potential biases will already have been ac- 

5 There is an interesting implication in blazar variability for the 
extragalactic gamma-ray background, which is discussed in jj4] 
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counted for in the way we have dealt with truncations 
in the Fioo,r plane. We have treated the truncation in 
tha t pla ne as empirical and accounted for it as discussed 
in 

In summary, most of these additional sources of er- 
ror are small and are important only for a small range of 
fluxes around the lowest fluxes, where the EP method has 
a larger uncertainty anyway, and for that reason these 
low fluxes are excluded from the fitting in our analysis 
(compare the end points of raw and corrected distribu- 
tions in Figures 3 and 4). 

4. TOTAL OUTPUT FROM BLAZARS AND THE 
EXTRAGALCTIC GAMMA-RAY BACKGROUND 

We can use the above results to calculate the total flux 
from blazars and the contribution of blazars to the EGB, 
defined here as the total extragalactic gamma-ray photon 
output The total output in gamma-ray photons from 
blazar sources with fluxes greater than a given Fioo, in 
terms of photons s _1 cm -2 sr -1 between 0.1 to 100 GeV, 
is 

POO 

J 7 (>F 100 )=/ ^'oo^i'ooXoo- (15) 

Integrating by parts the contribution to the EGB can be 
related directly to the cumulative distribution $(Fioo) 
which is the primary output of our procedure 

/>OC 

£y(> Fioo) = ^ioo $(Fioo) + / W 00 ) dF[ m . (16) 

Jf 100 

The advantage of using the latter equation is that it can 
give a step-by-step cumulative total contribution to the 
background instead of using analytic fits to the differ- 
ential or cumulative distributions obtained from binning 
the data. Figure shows X 1 {> Fioo) resulting from this 
integration down to flux Fioo = 5 x 10 -9 , with the to- 
tal output of the blazar population at fluxes probed by 
this analysis being X 7 (> Fioo = 5 x 10 _9 )=4.5 ± 0.5 
xlO -6 photons s _1 cm -2 sr -1 . Note that this includes 
the contribution from both detected blazars and those 
undetected above this flux owing to the truncation in 
the Fioo, T plane. Therefore, as expected it is more than 
the total contribution of blazars resolved by Fermi-LAT 
which is esti mated to be 4.1 ± 0.2 x 10 -6 photons s _1 
cm -2 sr~ 1 (jAbdo et all 2010c) PI 

In order to determine the contribution of blazars with 
Fioo < 5 x 10 -9 photons s _1 cm -2 sec -1 to the to- 
tal EGB we must extrapolate the flux distribution we 
have obtained to below this flux which cannot be unique 
and is more uncertain. We fit a power law to the faint 
end of <£>(Fioo) so that we can extend the integration of 
equation [16] to lower fluxes. Extending to zero flux we 
find that blazars in toto can produce a photon output 
of X 7 (> Fioo = 0)=8.5 (+6.3/-2.1) xl(T 6 photons s" 1 
cm -2 sr -1 . This large range is due to the uncertainty in 

6 This definition avoids the problem that individual instruments 
resolve a different fraction of sources, leading to different estimates 
for the fraction of the total extragalactic photon output that is 
unresolved. 

7 Actually the latter is what is attributed to point sources with 
test statistic value of TS > 25 which corresponds roughly to a 5cr 
detection. 



the faint end cumulative source counts slope, ultimately 
owing to the range of the correlation parameter /3, where 
the best fit value reported is for the middle of the 1<j 
faint end slope of $(Fioo)- 

This is to be compared with the total observed EGB 
of Xegb = 14.4 ± 1.9 x 10" 6 photons s" 1 cm" 2 sr -1 re- 
ported by Fermi^\ If the blazar population continues to 
have the fitted power law distribution to zero flux then 
it is clear that for our best fit parameters blazars can 
produce 59% of the observed EGB but this contribution 
could be as little as 39% or as much as all of the total 
extragalactic gamma-ray output of the Universe. This 
result is in agreement, albeit with a larger uncertainty, 
with the result in MA, where following the definition con- 
ventions here blazars extrapolated to zero flux are found 
to contribute 46% ± 10% of the EGB0 

It is, however, likely that blazars do not continue as a 
population with no change in the source counts slope to 
zero flux, since even a dim AGN of luminosity ~10 45 erg 
s _1 at redshift 3 would have a flux of ~10 -12 photons 
cm -2 sec -1 . If we only integrate equation[l6]to this lower 
flux limit, then we get the total blazar contribution to be 
£y(> Fioo = 1 x 10- 12 )=7.7 (+0.8/-1.2) xl(T 6 photons 
s -1 cm -2 sr -1 , which brings the upper limit estimate 
down to 66% of the total EGB. 

We can also obtain the energy intensity of the cumu- 
lative emission from blazars as 

rOO rOO 

^ 7 -blazars(> ^100 ) = / ^100 / dTE(F{ 00 ,T) G(F[ 00 ,T) 

(17) 

where for a simple power law spectrum we can relate the 
energy emitted between 0.1 and 100 GeV to the flux as 

e r — i i — io 3 ( 2-r ) 

^ - m * 100 x _ X l _ l03(1 _ r) MeV/photon, 

(18) 

except for T = 2 and T = 1 for which R(2) = In 10 3 /(1 - 
10" 3 ) - 6.9 and R(l) = (10 3 - 1)/Inl0 3 - 150 re- 
spectively. If we ignore the weak correlation between 
r and Fioo (set f3 = 0)_we get / 7 - b iazars(> ^ioo) = 
R x Xy(> .Fioo) where R is average value of R over 
the Gaussian distribution of T. We carry this average 
numerically and get the total (resolved and unresolved) 
/ 7 -biazars = 2.7(+3.1/ - 0.9) x 1(T 3 MeV cm" 2 sec -1 
sr -1 , integrating to zero flux taking into account the un- 
certainties above and the la uncertainty in the mean of 

We note that this analysis can not rule out blazars 

8 The Fermi collaboration papers divide this radiation into 
two parts, one from what is referred to as the contribution of 
resolved sources (the T sources = 4.1 ± 0.2 x 10 -6 photons s _1 
cm -2 sr -1 mentioned above), and a second "diffuse" component 
of Tegb- sources = 1-03 ± 0.17 x 10 -5 photons s _1 cm -2 sr -1 
Abdo et al. (2010d). However, the most relevant comparison is 
with the total of these two, because which sources are declared to 
be resolved is determined by a TS threshold, not a flux limit, and 
these are different due in part to the truncation in the i^ioo, T plane 
and the varying Galactic diffuse level. 

9 The total point source diffuse emission and EGB intensity pre- 
sented in Table 6 of MA have the contribution of resolved Fermi- 
LAT sources removed, so for direct comparison to the results pre- 
sented here the total Fermi-LAT resolved source contribution of 
4.1 ± 0.2 x 10 -6 photons s _1 cm -2 sr -1 must be added to both 
before the ratio is taken. 
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as the sole significant contributor to the EGB, although 
the best fit value does not favor this bei ng the case. 
The s pectral index of the EGB of ~2.4 ([Abdo et all 
l2010df ) is consistent with the mean photon index of 
the blaza rs as determined here and in MA. In a sim- 
ilar vein, IVenters fc Pavildoul (|2011[ ) have shown that 
the spectrum of the EGB is consistent with a blazar 
origin. Several auth ors (e.g. iStecker fc Ventersl 1201 lb 
lAbazajian et al.ll2010f ) have suggested that blazars could 
be the primary sourc e of the EGB, wh il e the results 
presented in MA and iMalvshev fc Hoggl (j2011f ) would 
favor the primary source being something else. Other 
possible source populations include starforming galax- 
ies, which have been recognize d as a possible ma j or con - 
tributor to the EG B by e.g. IStecker fc Ventersl ([201 



iFields et all ([2010D . and lLacki et all (120111). a lthough 
this has bee n disputed b y iMakiva et aL I (j20ll . radio 
galax ies fe.g llnoue 11201 ih, and other non-blazar AGN 

(e.g. llnoue fc Tot anil I2009L l20TTh. 

According to IStecker fc Salamonl (j 19961 ) , to the extent 
that faint blazars are more likely to be observed by in- 
struments such as the Fermi-L AT if they are in the flaring 
state rather than the quiescent state, then the observed 
blazars should have a different mean mean photon in- 
dex than the EGB, were the EGB to be made primarily 
from quiescent state blazars, under the assumption that 
blazars in the flaring state have a different spectrum than 
in the quiescent state. As the reconstructed mean pho- 
ton index here of the Fermi-LAT observed population is 
close to that of the EGB, and there is only a weak rela- 
tion and correlation between flux and photon index, this 
would imply that at least one of the following must be the 
case: a) there is not a significant bias in the Fermi-LAT 
toward detecting blazars in the flaring state, b) quiescent 
blazars do not form the bulk of the EGB, or c) flaring 
and quiescent blazars have, en masse, roughly the same 
photon index distributions. 

5. DISCUSSION 

We have used a rigorous method to calculate the intrin- 
sic distributions in flux (known commonly as the source 
counts or the LogA^-Log6' relation) and photon index 
of Fermi-LAT blazars directly from the observed ones 
without any assumptions or reliance on extensive sim- 
ulations. This method features a robust accounting for 
the pronounced data truncation introduced by the selec- 
tion biases inherent in the observations, and addresses 
the possible correlation between the variables. The ac- 
curacy of the methods used here are demonstrated in the 
Appendix using a simulated data set with known distri- 
butions. A summary of the best fit correlations between 
photon index and flux, and the best fit parameters de- 
scribing the inherent distributions of flux and photon in- 
dex, of the observed data are presented in Table 1 along 
with the values obtained by MA. We have obtained the 
distributions of flux and photon index of blazars consid- 
ering the major data truncation arising from Fermi-LAT 
observations. More subtle issues affecting the distribu- 
tions we have derived, especially the photon index dis- 
tribution, may arise due to the finite bandwidth of the 
Fermi-LAT and lack of complete knowledge of the ob- 
jects' spectra over a large energy range and deviations 
from simple power laws. However the Fermi-LAT band- 
width is sufficiently large that the contribution of sources 




F 100 [photons cm ; 



Fig. 7. — Estimate of the cumulative number of photons between 
0.1 and 100 GeV above a given Fioo from blazars, X 7 _bi a zars(> 
^100 ) from equation 1161 shown with error bars resulting from the 
lcr range of the correlation parameter f3. The bottom dotted hor- 
izontal line shows the level of the EGB as measured by Fermi 
( Abdo et al. 2010d), with Fermi resolved point sources removed. 
The top dashed horizontal line shows the EGB, ie the total ex- 
tragalactic gamma-ray output (Xegb) as denned here. The best 
fit total con trib ution of blazars to Xegb, obtained by integrat- 
ing equation 1161 to zero flux, is 59%, but our analysis cannot rule 
out blazars, integrated to arbitrarily low flux, forming as little as 
33% or as much as all of the total extragalactic gamma-ray out- 
put, due to the large range of uncertainty in the faint end source 
counts slope, ultimately owing to the uncertainty in the correlation 
parameter (3. 

which peak outside of this range to the source counts and 
the EGB in this energy range will be small. 

We find that the photon index and flux show a slight 
correlation, although this correlation is of marginal sig- 
nificance. This indicates that the intrinsic luminosities 
and photon indexes are correlated only weakly. The com- 
parison of the intrinsic and raw observed distributions 
show clearly the substantial effects of the observational 
bias. The intrinsic differential counts can be fitted ade- 
quately by a broken power law and the photon index ap- 
pears to have a intrinsic Gaussian distribution. We also 
find that in general the values reported here are consis- 
tent with those reported in MA for the power law slopes 
of the flux distribution V{^ioo) and the distributions of 
photon index ft(r), although the allowed range of the cor- 
relation parameter /? here allows for wider uncertainty in 
these values in some cases. We do note a discrepancy at 
the la level for the faint end slope of the FSRQ source 
counts. 

Using the bias free distributions we calculated the total 
cumulative contribution of blazars to the EGB as a func- 
tion of flux. We obtain this directly from the cumulative 
flux distribution which is the main output of the meth- 
ods used. Under the assumption that the distribution 
of blazars continues to arbitrarily low flux, we find the 
best fit contribution of blazars to the total extragalactic 
gamma-ray radiation in the range from 0.1 to 100 GeV to 
be at the level of 59%, although this analysis cannot rule 
out blazars producing as little as 39% or as much as all 
of the total extragalactic gamma-ray output. This result 
is in agreement with the result in MA, although with a 
larger uncertainty. The significant uncertainties reported 
here for the source count slopes and the contribution of 
blazars to the EGB are ultimately due to the allowed 
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TABLE 1 

Best fit parameters for Fermi-LAT blazar intrinsic distributions, as calculated in this work and in MA (|Abdo et al.II201"0bT) 





n 


(3 a 


Wlabove 


Fbreak 


Wlbelow 


M e 


a { 


Blazars s (this work) 
BlazarsS (MA) 


352 
352 


0.02±0.08 


-2.37±0.13 
-2.48±0.13 


7.0±0.2 
7.39±1.01 


-1.70±0.26 
-1.57±0.09 


2.41±0.13 
2.37±0.02 


0.25±0.03 
0.28±0.01 


BL Lacs (this work) 
BL Lacs (MA) 


163 
163 


0.04±0.09 


-2.55±0.17 
-2.74±0.30 


6.5 ±0.5 
6.77±1.30 


-1.61±0.27 
-1.72±0.14 


2.13±0.13 
2.18±0.02 


0.24±0.02 
0.23±0.01 


FSRQs (this work) 
FSRQs (MA) 


161 
161 


-0.11±0.06 


-2.22±0.09 
-2.41±0.16 


5.1±2.0 
6.12±1.30 


-1.62±0.46 
-0.70±0.30 


2.52±0.08 
2.48±0.02 


0.17±0.02 
0.18±0.01 



a The correlation between photon index F and Log flux Fioo- See Equation [2] and fl3.ll A higher value of f3 (i.e. more positive correlation between 
flux and photon index absolute value) moves the mean of the photon index distribution down to lower photon index absolute value (lower /x) and 
makes the faint end source counts slope less steep (less negative rribeiow), while a lower value of (3 has the opposite effect. All values reported for 
this work include the full range of results and their uncertainties when considering the lcr range of (3. 
b The power law of the intrinsic flux distribution t/^-Fioo) at fluxes above the break in the distribution. See Equation [9] 

c The flux at which the power law break in ^-^ioo) occurs, in units of 10 -8 photons cm -2 sec -1 . We present the value even though the precise 

location of the break is not important for the analysis in this work. The value of Fb rea k can be obtained by equations 1101 and II II 

d The power law of the intrinsic flux distribution ^Fioo) at fluxes below the break. See Equation [9] 

e The mean of the Gaussian fit to the intrinsic photon index distribution h{T). 

f The lcr width of the Gaussian fit to the intrinsic photon index distribution h{F). 

g Including all FQRQs, BL Lacs, and 28 of unidentified type. 



range of the correlation parameter j3. As discussed in 
the Appendix, the method applied to the (uncorrelated) 
simulated data also manifests a significant uncertainty on 
/?, that translates into the corresponding uncertainty for 
the faint-end slope of the source count distribution. This 
is important as it ultimately governs the contribution of 
blazars to the EGB. We note that a similar scenario (i.e. 
absence of a correlation between photon index and flux) 
might characterize the real data in view of the results 
reported in the previous sections and their similarity to 
the results obtained using simulated uncorrelated data. 
As shown in the Appendix, larger samples are required 
to narrow down the uncertainty on the correlation pa- 
rameter /3. 

If, as could be expected, the flux distribution flattens 
at fluxes below ~10 -12 photons cm -2 sec -1 , the inte- 



grated contribution will be significantly lower than for a 
naive extrapolation to zero flux. This is also modulo any 
change in the power law slope of the source counts below 
the fluxes probed in this analysis, which might arise due 
to luminosity and/or density evolution with redshift. A 
full accounting for the possible evolution in the blazar 
population using a sample with redshift determinations 
will be presented in a forthcoming work. 
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APPENDIX 

Here we test the methods used in this paper with a simulated Monte Carlo data set provided by the Fermi-LAT 
collaboration. The data set is a realization of a set of simulations discussed in MA. For each Monte Carlo realization 
20000 sources were placed isotropically on the sky according to assumed distributions of flux (broken power law) and 
photon index (Gaussian). Instrumental and observational effects on detection were applied to this data, resulting in a 
catalog of 486 sources with detection TS of at least 50, where TS is the test statistic as discussed in £[2j 

Figure [8] shows the fluxes and photon indexes for the simulated data set, along with a curve approximating the 
observation truncation in the Fi o — T plane, which we take as the limiting flux for any object at a given photon index, 
and the limiting photon index for any object at a given flux, for inclusion in the relevant associated sets. Since there 
is some uncertainty in the detection threshold values of fluxes and indexes we carry our analysis first assuming the 
individual limiting fluxes for each source to be Fu m = Fioo / and then by using a simple curve to define the 
truncation boundary as shown in Figure [8] This is a more conservative assumption and few sources are excluded but 
it insures the completeness of the sample. We experiment with moving the curve to the right and down (eliminating 
more sources in the edges of the sample) until we do not notice any change in the result. As described below this 
reproduces the input data accurately. We carry out the same procedure for the real data. 

While the raw data obtained from these simulations show a strong correlation between flux and index (ignoring the 
truncation we obtain the correlation parameter /3 raw = 0.53 ±0.03 defined in equation [2]) , our method shows that once 
the effects of truncation are accounted for the correlation disappears and we get a correlation parameter j3 consistent 
with zero, in agreement with that of the input data. Figure [9] shows the values of j3 vs r for both the raw data and 
with the truncation accounted for. 

Figure [10] shows the observed and reconstructed intrinsic differential distribution of flux ^Fioo) for the simulated data 
set, along with the known intrinsic distribution. Figure fTTI shows the observed and reconstructed intrinsic differential 
distribution of photon index HY) for the simulated data set, along with the known intrinsic distribution. Table 2 
shows the full range of values for the reconstructed intrinsic distributions for the simulated data, including the effects 
of the entire la range of the correlation parameter /3, along with the known input distributions. It is seen that these 
methods successfully reproduce the input intrinsic distributions from a highly truncated data set. 

Here we can also address the question of whether errors in the cumulative distribution propagate in a significantly 
correlated way to the differential distributions. First we take random samples of size n from the simulated data set and 
add 0.3 times the flux of each object to itself in half of each of the samples and subtract 0.3 times the flux of each object 
from itself in the other half. This factor of 0.3 reflects the largest typical reported uncertainties in the Fermi-LAT flux 
measurements that we use. The effect of doing so on the cumulative flux distribution is shown in Figure [12] while the 
effect propagated through to the cumulative total number of photons (i.e. the contribution to the EGB), determined in 
the manner of equation [16] is shown in Figure [13] In both figures the solid curve shows the cumulative flux distribution 
<I>(Fioo) for n=0, the base case of no changes. The dashed curves show n=10 and the dotted curves are for n=100, the 
later representing almost one quarter of the objects in the sample. In these cases the differences in the cumulative flux 
distribution and the cumulative number of photons from the base case are negligible. Then, for a more realistic but 
perhaps extreme case, we alter the flux of all of the objects with alterations distributed such that those objects with 
the lowest fluxes have their fluxes altered by the highest typical reported measurement errors of 0.3 times the flux, 
while those with higher fluxes have lower errors, with the alteration proportional to the ratio of the difference of the 
logarithm of the flux and that of the maximum flux in the sample, with positive alterations for half the objects and 
negative alterations for half. The resulting cumulative flux distribution and cumulative number of photons is shown 
by the dash-dot curves in Figures [12] and [13] Even then the change to the cumulative flux distribution is small and 
the added uncertainty introduced into the fitted differential distribution and the cumulative total number of photons 
if this case is considered relative to the base case is small compared to the uncertainty resulting from considering the 
extremal values of the correlation parameter f3 and other sources of uncertainty considered in this work (compare with 
Figures [3] and [7]). 

We have also examined the effect of increasing the data set size on the uncertainty range determined for the 
correlation parameter j3 for the method employed in this work. A second simulated data set provided by the Fermi- 
LAT collaboration consisting of a catalog of ~6 times as many (3018) objects resulted in a la range for j3 that was 
approximately half as large (/3=-0.01=b -0.04) as with the 486 object set. As the uncertainty range in the value of /3 is 
the major driver in the total uncertainty of the fitted distribution parameters, we can expect a significant reduction 
in uncertainty levels for the distribution parameters of the real data with a future five year Fermi-LAT extragalactic 
catalog consisting of ~1500 blazars as opposed to the 352 here. 
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Fig. 8. — Flux and photon index for the 486 sources in the simulated Fermi-LAT data set, along with the curve (solid line) used to 
specify the observation truncation in the Fioo,T plane. As with the real blazar data, moving the cutoff to the left dashed line has a large 
effect on the results, but moving it to the right dashed line has a negligible effect. 
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Fig. 9. — Correlation factor (3 versus test statistic r for a photon index and flux correlation of the form given in Equation [2] for the 486 
sources in the simulated Fermi-LAT data set. The solid curve show the results for the method employed here with the cutoff curve shown 
in Figure [8] while the dotted curve shows the results for the raw data. 



TABLE 2 

Simulated blazar intrinsic distributions, as calculated in this work and as known from the inputs to the simulated dataset 

oa ryi b 771 c ryi d ..e _f 
P ru above r break ru below A* v 



Determined distributions 8 0.02±0.07 -2.41±0.11 8.0±1.0 -1.61±0.27 2.38±0.9 0.34±0.05 

Known input distributions^ -2.49 6.6 -1.59 2.37 0.28 



a The correlation between photon index F and Log flux Fioo- See Equation [2] and fc|3.1| 

b The power law of the intrinsic flux distribution ^-^ioo) at fluxes above the break in the distribution. See Euqation[9] 
c The flux at which the power law break in ^-Fioo) occurs, in units of 10 — 8 photons cm -2 sec -1 . 
d The power law of the intrinsic flux distribution ^Fioo) at fluxes below the break. See Equation [9] 

e The mean of the Gaussian fit to the intrinsic photon index distribution h(T). For the analysis here this includes the full range of results and their 
uncertainties when considering the la range of (3. 

f The la width of the Gaussian fit to the intrinsic photon index distribution h(T). 

g The simulated observed blazar data set is provided by the Fermi collaboration, and has 486 objects. All values for the determined distributions 
reported here include the full range of results and their uncertainties when considering the la range of (3. 
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Fig. 10. — Observed (diamonds) and reconstructed ft = intrinsic (stars) differential distribution of flux t/^-^iOo) for the 486 sources in 
the simulated Fermi-LAT data set. The data have intrinsic power law slope distributions of -2.49 and -1.59 above and below the break, 
respectively, which are plotted. The normalization of V^ioo) here is arbitrary. 




photon index 



Fig. 11. — Observed (diamonds) and reconstructed (3 = intrinsic (stars) differential distribution of photon index h(T) for the 486 
sources in the simulated Fermi-LAT data set. The data have an intrinsic Gaussian distribution with a mean of 2.37 and la width of 0.28, 
which is shown by the dashed curve. The solid curve is the best fit Gaussian to the stars, which differs only slightly from the dashed curve. 
The normalization of h(T) here is arbitrary. 
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Fig. 12. — Reconstructed (3 = intrinsic cumulative distribution of flux ^Fioo) for the 486 sources in the simulated Fermi-LAT data set. 
The normalization of $(Fioo) here is arbitrary. The solid curve shows the fluxes as simulated, while the dashed (n=10) and dotted (n=100) 
curves show the results if a number n of those fluxes are altered in such a way that half of the altered fluxes are increased by 30% and 
half are decreased by 30%, values representing the largest typical reported uncertainties for the flux measurements used in this analysis. 
The dash-dot curve shows the result for a more realistic case where the fluxes of all of the objects are altered in the manner described in 
the text. As evident in all cases the added uncertainty introduced in the cumulative and fitted differential distribution is small compared 
to the uncertainty resulting from considering the extremal values of the correlation parameter (3 or other uncertainties considered in this 
analysis. 
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Fig. 13. — Same as Figure \12\ except the reconstructed f$ = estimate of the cumulative total number of photons between 0.1 and 
100 GeV above a given Fioo from blazars, X 7 _bi a zars(> ^ioo), from equation 1161 The normalization of X 7 _bi a zars(> ^ioo) here for the 
simulated data set is arbitrary. 



